/****************************************************************************
 * Copyright (C) 2009-2010 SciTouch LLC
 * 
 * This file is part of Indigo toolkit.
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 3 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 ***************************************************************************/

#include "molecule/molecule_tautomer.h"

#include "base_cpp/queue.h"
#include "molecule/molecule.h"
#include "molecule/molecule_substructure_matcher.h"

// Tautomer superstrucure - molecule with all bonds that can appear in tautomer
// This structure is needed to enumerate all submolecule of a tautomer

// Tautomer template rules:
// A - any atom
// G - heteroatom: N, O, S, P, As, Sb, Se, Te
// Tautomer bond appear/disappear template:
//    G=C   G-H        <-->  G-C-G-H
//    (A=A-)nG=C  G-H  <-->  H-(A-A=)nG-C-G
// Atoms can't have any other double/triple bonds
// Lets say thet C can emit bond and G with H can accept bond

Molecule& TautomerSuperStructure::getSuperMolecule (Array<int> **mapping, Array<int> **inv_mapping)
{
   if (mapping != 0)
      *mapping = &_mapping;

   if (inv_mapping != 0)
      *inv_mapping = &_inv_mapping;

   return _superMolecule;
}

int TautomerSuperStructure::getSubgraphType (const Array<int> &vertices, const Array<int> &edges)
{
   // For any atoms number of attached bonds must be 0 or 1

   QS_DEF(Array<int>, per_vertex_attached_bonds);
   per_vertex_attached_bonds.clear_resize(_superMolecule.vertexEnd());
   for (int i = 0; i < vertices.size(); i++)
   {
      int v = vertices[i];
      per_vertex_attached_bonds[v] = 0;
   }

   int attached_bonds = 0;
   for (int i = 0; i < edges.size(); i++)
   {
      int edge_index = edges[i];

      if (!_isBondAttachedArray[edge_index])
         continue;

      const Edge &edge = _superMolecule.getEdge(edge_index);
      per_vertex_attached_bonds[edge.beg]++;
      per_vertex_attached_bonds[edge.end]++;
      if (per_vertex_attached_bonds[edge.beg] > 1 || per_vertex_attached_bonds[edge.end] > 1)
         return NONE;
      attached_bonds++;
   }

   if (attached_bonds == 0)
      return ORIGINAL;
   return TAUTOMER;
}

TautomerSuperStructure::TautomerSuperStructure () : 
   TL_CP_GET(_atomsEmitBond),
   TL_CP_GET(_atomsAcceptBond),
   TL_CP_GET(_isBondAttachedArray),
   TL_CP_GET(_superMolecule),
   TL_CP_GET(_mapping),
   TL_CP_GET(_inv_mapping)
{
}

void TautomerSuperStructure::setBaseMolecule (const Molecule &baseMolecule)
{
   _superMolecule.cloneMolecule(baseMolecule, &_mapping, &_inv_mapping);
}

void TautomerSuperStructure::build ()
{
   // Collect atom properties
   _collectAtomProperties();

   // Detect distances from _atomsEmitBond elements to _atomsAcceptBond elements
   QS_DEF(Array<int>, distancesMatrix);
   distancesMatrix.resize(_atomsEmitBond.size() * _atomsAcceptBond.size());
   for (int i = 0; i < _atomsEmitBond.size(); i++)
   {
      int *result = distancesMatrix.ptr() + _atomsAcceptBond.size() * i;
      _findMinDistance(_atomsEmitBond[i], 6, _atomsAcceptBond, result);
   }


   QS_DEF(Array<int>, attachedBonds);
   attachedBonds.clear();
   for (int i = 0; i < _atomsEmitBond.size(); i++)
      for (int j = 0; j < _atomsAcceptBond.size(); j++)
      {
         int v1 = _atomsEmitBond[i];
         int v2 = _atomsAcceptBond[j];
         if (_superMolecule.findEdgeIndex(v1, v2) != -1)
            continue;
         // Check new loop size: 5 or 6
         int size = distancesMatrix[_atomsAcceptBond.size() * i + j];
         if (size != 4 && size != 5)
            continue;

         int bond = _superMolecule.addBond(v1, v2, BOND_ZEROED);
         attachedBonds.push(bond);
      }

   _isBondAttachedArray.resize(_superMolecule.edgeEnd());
   _isBondAttachedArray.zerofill();
   for (int i = 0; i < attachedBonds.size(); i++)
      _isBondAttachedArray[attachedBonds[i]] = true;
}

void TautomerSuperStructure::restore ()
{
   for (int i = 0; i < _isBondAttachedArray.size(); i++)
      if (_isBondAttachedArray[i])
         _superMolecule.removeEdge(i);
}

bool TautomerSuperStructure::_isEmitHeteroatom (const Atom &atom)
{
   static const int heteroatoms[] = { ELEM_O, ELEM_N };

   for (unsigned int i = 0; i < NELEM(heteroatoms); i++)
      if (atom.label == heteroatoms[i])
         return true;
   return false;
}

void TautomerSuperStructure::_getDoubleBondsCount (int i, int &double_count, 
                                                  int &arom_count)
{
   double_count = 0;
   arom_count = 0;

   const Vertex &vertex = _superMolecule.getVertex(i);
   for (int j = vertex.neiBegin(); j != vertex.neiEnd(); j = vertex.neiNext(j))
   {
      const Bond &bond = _superMolecule.getBond(vertex.neiEdge(j));
      if (bond.type == BOND_DOUBLE || bond.type == BOND_TRIPLE)
         double_count++;
      else if (bond.type == BOND_AROMATIC)
         arom_count++;
   }
}

void TautomerSuperStructure::_collectAtomProperties (void)
{
   QS_DEF(Array<int>, h_counters);
   MoleculeSubstructureMatcher::calcTargetHCounters(_superMolecule, h_counters);

   _atomsAcceptBond.clear();
   _atomsEmitBond.clear();
   for (int i = _superMolecule.vertexBegin(); 
        i != _superMolecule.vertexEnd(); 
        i = _superMolecule.vertexNext(i))
   {
      const Atom &atom = _superMolecule.getAtom(i);
      bool have_hydrogen = (h_counters[i] != 0);
      int  double_bonds_count, arom_bonds_count;

      _getDoubleBondsCount(i, double_bonds_count, arom_bonds_count);
      if (arom_bonds_count != 0)
         have_hydrogen = true;

      // Detect atom type
      if (double_bonds_count == 0 && have_hydrogen && Element::isHeteroatom(atom.label))
         _atomsAcceptBond.push(i);

      // Aromatic bond can be double
      if ((double_bonds_count == 1 || arom_bonds_count != 0) && atom.label == ELEM_C)
      {
         int nei_heteroatoms = _hetroatomsCount(i);
         if (nei_heteroatoms != 1)
            continue;

         // Check all possible double bond neighbors to have emit atom type
         const Vertex &vertex = _superMolecule.getVertex(i);
         for (int j = vertex.neiBegin(); j != vertex.neiEnd(); j = vertex.neiNext(j))
         {
            // Double bond can jump from non-heteroatom to heteroatom 
            // that why only existance of heteroatom and existance of any attached double bond
            // is required (double bond count already has been checked above)

            const Atom &neighbourAtom = _superMolecule.getAtom(vertex.neiVertex(j));
            if (!_isEmitHeteroatom(neighbourAtom))
               continue;

            int nei_double_bonds_count, nei_arom_bonds_count;
            _getDoubleBondsCount(vertex.neiVertex(j), nei_double_bonds_count, nei_arom_bonds_count);
            if (nei_double_bonds_count <= 1)
            {
               _atomsEmitBond.push(i);
               break;
            }
         }
      }
   }
}

// Find minimum distance between source vertex and vertices from dest array
// and check tautomer chain property:
//    * only one bond can be in ring
//    * if no one bond is in ring then only one bond can be double and no one can be triple
void TautomerSuperStructure::_findMinDistance (int source, int maxDist, Array<int> &dest, int *result)
{
   QS_DEF(Array<int>, distances);
   QS_DEF(Array<int>, parents);
   distances.resize(_superMolecule.vertexEnd());
   parents.resize(_superMolecule.vertexEnd());

   // Fill distances by infinity
   const int INFINITY = 100000;
   for (int j = 0; j < distances.size(); j++)
      distances[j] = INFINITY;
   QS_DEF(Queue<int>, front);
   front.clear();
   front.setLength(_superMolecule.vertexEnd());

   distances[source] = 0;
   parents[source] = -1;
   front.push(source);
   while (!front.isEmpty())
   {
      int active = front.pop();
      if (distances[active] == maxDist)
         break;

      const Vertex &vertex = _superMolecule.getVertex(active);
      for (int j = vertex.neiBegin(); j != vertex.neiEnd(); j = vertex.neiNext(j))
      {
         int vn = vertex.neiVertex(j);
         if (distances[vn] == INFINITY)
         {
            distances[vn] = distances[active] + 1;
            parents[vn] = active;
            front.push(vn);
         }
      }
   }
   for (int j = 0; j < dest.size(); j++)
   {
      // Check chain
      if (distances[dest[j]] != INFINITY)
      {
         int inRingCount = 0;
         int doubleBondsCount = 0, tripleBondsCount = 0;

         int cur = dest[j];
         int prev = parents[cur];
         while (prev != -1)
         {
            int edge_idx = _superMolecule.findEdgeIndex(cur, prev);
            if (_superMolecule.bondInRing(edge_idx))   
               inRingCount++;
            if (_superMolecule.getBond(edge_idx).type == BOND_DOUBLE)   
               doubleBondsCount++;
            if (_superMolecule.getBond(edge_idx).type == BOND_TRIPLE)   
               tripleBondsCount++;

            cur = prev;
            prev = parents[prev];
         }

         if (inRingCount > 1)
            distances[dest[j]] = 2 * INFINITY;
         else if (inRingCount == 0)
            if (doubleBondsCount > 1 || tripleBondsCount > 0)
               distances[dest[j]] = 2 * INFINITY;
      }
      result[j] = distances[dest[j]];
   }
}

int TautomerSuperStructure::_hetroatomsCount (int idx)
{
   const Vertex &vertex = _superMolecule.getVertex(idx);
   int count = 0;
   for (int i = vertex.neiBegin(); i != vertex.neiEnd(); i = vertex.neiNext(i))
   {
      const Atom &atom = _superMolecule.getAtom(vertex.neiVertex(i));
      if (Element::isHeteroatom(atom.label))
         count++;
   }
   return count;
}
